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Abstract. We revisit the totally asymmetric simple exclusion process with open 
boundaries (TASEP), focussing on the recent discovery by de Gier and Essler that the 
model has a dynamical transition along a nontrivial line in the phase diagram. This 
line coincides neither with any change in the steady-state properties of the TASEP, 
nor the corresponding line predicted by domain wall theory. We provide numerical 
evidence that the TASEP indeed has a dynamical transition along the de Gier-Essler 
line, finding that the most convincing evidence was obtained from Density Matrix 
Renormalisation Group (DMRG) calculations. By contrast, we find that the dynamical 
transition is rather hard to see in direct Monte Carlo simulations of the TASEP. We 
furthermore discuss in general terms scenarios that admit a distinction between static 
and dynamic phase behaviour. 
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1. Introduction 

As a rare example of a class of exactly-solvable, nonequilibrium interacting particle 
systems, asymmetric simple exclusion processes of various types have found favour with 
researchers in statistical mechanics, mathematical physics and probability theory [IHI] . 
These systems comprise hard-core particles hopping in a preferred direction on a one- 
dimensional lattice and have been used to describe systems as diverse as traffic flow [5] , 
the dynamics of ribosomes [6] and molecular motors [318], the flow of hydrocarbons 
through a zeolite pore [9], the growth of a fungal filament [10J, and in queueing 
theory [HJ[12]. The physical interest lies mainly in the rich phase behaviour that arises 
as a direct consequence of being driven away from equilibrium. 

Of particular interest are the open boundary cases. Here the system can be thought 
of as being placed between two boundary reservoirs, generally of different densities. The 
two reservoirs enforce a steady current across the lattice, and therewith a nonequilibrium 
steady state. It was first argued by Krug [13] that as the boundary densities are varied 
nonequilibrium phase transitions may occur in which steady state bulk quantities — 
such as the mean current or density — exhibit nonanalyticities. Such phase transitions 
were seen explicitly in first a mean-field approximation [HUH] and then in the exact 
solution [151 US] of the totally asymmetric exclusion process with open boundaries, 
hereafter abbreviated as TASEP, and of related processes p0|. Since then arguments 
have been developed to predict the phase diagram of more general one-dimensional 
driven diffusive systems [T7I4T9] . 

Our interest in this work is in the distinction between two phase diagrams for 
the TASEP: one of which characterises the steady-state behaviour, and the other the 



dynamics. The static phase diagram [T5j[T6] is shown in Fig. 1(a) where the left boundary 
density is a and the right is 1 — (3. There are three possible phases in the steady state: 
a low density (LD) phase where the bulk density is controlled by the left boundary and 
is equal to a; a high density (HD) phase where the bulk density is controlled by the 
right boundary and is equal to 1 — /3; a maximal current (MC) phase where the bulk 
density is |. The high and low density phases are further divided into subphases (e.g. 
LD1, LDII) according to the functional form of the spatial decay of the density profile 
to the bulk value near the non-controlling boundary. In these subphases the lengthscale 
over which the decay is observable remains finite as the system size is increased. This 
change at the subphase boundaries in the form of the decay is thus not a true phase 
transition in the thermodynamic sense. 

More recently, de Gier and Essler have performed an exact analysis of the ASEP's 
dynamics [20H22] and derived the the longest relaxation times of the system by 
calculating the gap in the spectrum using the Bethe ansatz. The analysis builds on 
directly related work on the XXZ spin chain with nondiagonal boundary fields [23,;24j. 
The dynamical phase diagram they obtain, in which phases are associated with different 



functional forms of the longest relaxation time, is illustrated in figure 1(b) The same 



phases (high density, low density and maximal current) are found as in the static phase 
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Figure 1. (a) Static and |(b)| dynamic phase diagrams for the TASEP. Solid lines 



indicate thermodynamic phase transitions at which the current and bulk density are 
nonanalytic. Dotted lines indicate subphase boundaries. In the static case, the density 
profile a finite distance from the boundary is nonanalytic. In the dynamic case, the 
longest relaxation time exhibits a nonanalyticity. 



diagram, however, a hitherto unexpected dynamical transition line which subdivides the 
low density and high density phases is now apparent. This line which we shall refer to 
as the de Gier-Essler (dGE) line replaces the subphase boundaries ata = l/2,/3<l/2 
(3 = 1/2, a < 1/2 in the static phase diagram and, for example, subdivides the low 
density region into LDI', LDII'. 

Although there is nothing to rule out the prediction of dGE of a dynamical 
transition at a distinct location to any static transition, the result came as something of 
a surprise. This is because the static phase diagram had been successfully interpreted in 
terms of an effective, dynamical theory thought to be relevant for late-time dynamics, 
referred to as domain wall theory (DWT). We shall review DWT more fully in 
Section 12.41 For the moment we note that DWT correctly predicts the static subphase 
boundary and the associated change in the decay of the density profile, but does not 
predict the dGE line. Therefore, it was previously thought that a dynamical transition 
also occurred at the subphase boundary and initial calculations of relaxation times 
appeared to be consistent with this [25||26]. 

Our primary goal in this work is to establish numerically that a dynamical transition 
does indeed occur along the dGE line rather than at the static subphase boundary. We 
used two distinct approaches. First, we carried out direct Monte Carlo simulations of 
the TASEP dynamics and identified the dominant transient in three time-dependent 
observables. We find that these simulations are consistent with a dynamical transition 
at the dGE line but are not sufficiently accurate to rule out the scenario of the dynamical 
transition occurring at the subphase boundary. 
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We instead turn to a Density Matrix Renormalisation Group (DMRG) approach [27J 
to acquire convincing evidence that a dynamical transition occurs at the dGE line. The 
DMRG is an approximate means to obtain the lowest-lying energies and eigenstates of 
model Hamiltonians with short-range interactions. Although originally developed in the 
context of the Hubbard and Heisenberg models [28] (see also [29] for a comprehensive 
review), it has also been applied to nonequilibrium processes such as reaction-diffusion 
models [30J. In particular, it has also been used to investigate the spectrum of the 
TASEP [26J. However, the dynamical transition was not evident from the data presented 
in that work — perhaps because it had not been predicted at that time. By revisiting 
this approach, we obtain estimates of the longest relaxation time that allow us to rule 
out the domain wall theory prediction for the dynamical transition. 

The paper is organised as follows. We begin in Section [2] by recalling the definition 
of the TASEP, and by reviewing in more detail the static and dynamic phase diagrams 
discussed above. Then, in Section [3j we present our Monte Carlo simulation data, 
followed by the DMRG results in Section HI In Section [5] we return to more general 
questions about the distinction between the static and dynamic phase diagram with 
particular reference to different theoretical approaches to the TASEP. We conclude in 
Section [6] with some open questions. 

2. Model definition and phase diagrams 

Although we have alluded to the basic properties of the totally asymmetric simple 
exclusion process (TASEP) in the introduction, in the interest of a self-contained 
presentation we provide in this section a precise model definition and full details of 
the static and dynamic phase diagrams. 

2.1. Model definition 

The TASEP describes the biased diffusion of particles on a one-dimensional lattice with 
L sites. No more than one particle can occupy a given site, and overtaking is prohibited. 
The stochastic dynamics are expressed in terms of transition rates, for example a particle 
on a site i in the bulk hops to the right as a Poisson process with rate 1, but only if site 
i + 1 is unoccupied: At the left boundary only particle influx takes place, with rate a, 
and at the right boundary only particle outflux takes place, with rate (3, as shown in 
Fig. [2j The corresponding reservoir densities are a at the left and 1 — /3 at the right. 

2.2. Static phase diagram 

As noted in the introduction, the TASEP static phase diagram can be determined 
from, for example, exact expressions for the current and density profile in the steady 
state [T5|IT6]. The current takes three functional forms according to whether it is limited 
by a slow insertion rate (LD: a < f3, a < |), by a slow exit rate (HD: (3 < a, (3 < |) 
or by the exclusion interaction in the bulk (MC: a > > \). This behaviour of 
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Figure 2. The dynamics of the totally asymmetric simple exclusion process (TASEP) 
with open boundaries. Arrows show moves that may take place, and the labels indicate 
the corresponding stochastic rates. 



the current leads to the identification of the thermodynamic phases separated by solid 
lines in Fig. 1(a) Along these lines there are nonanalyticities in both the current and 
the density far from either boundary. The relevant expressions for these quantities are 



shown on Fig. l(a^ 



The density near one of the boundaries is nonanalytic along the lines a = ~ and 
f3 — \ in the HD and LD phases respectively. These subphases are shown dotted in 



Fig. 1(a) Inspection of the functional form of the density profile reveals that these 
are not true thermodynamic phase transitions, in the sense that the deviation from 
the bulk value extends only a finite distance into the bulk, and thus contributes only 
subextensively to the nonequilibrium analogue of a free energy (see e.g., [3T] for the 
definition of such a quantity). To be more explicit, consider for example the behaviour 
near the right boundary in the low-density phase. Here the bulk density is p = a. 
The mean occupancy of the lattice site positioned a distance j from the right boundary 
approaches in the thermodynamic limit L — > oo the form [HfT5|[T6] 

. U + Cl (/3)(^y «</3<i (LDI) 

3 \ a + cnfrfl&ffit- \<fl (LDII) 

In these expressions, cj and en are functions of the boundary parameters that we leave 
unspecified here so as to focus on the lengthscale of the exponential decay from the 
right boundary. As /3 is increased from zero, the decay length increases until = |. 
Then the decay length is constant, and the exponential is modulated by a power law. 
The density profile at the left boundary exhibits the same kind of nonanalyticity in the 
high-density phase as a is increased through ~ as a consequence of the particle-hole 
symmetry, p^i(a, 0) = 1 — pL-i((3,a), exhibited by the model. 



2.3. Dynamical phase diagram 

The dynamic phase diagram is obtained by examining a quantity referred to as the 
gap by de Gier and Essler J2DH22]- This is simply the largest nonzero eigenvalue of 
the transition matrix governing the continuous-time stochastic dynamics of the TASEP. 
More formally, one starts with the master equation 

±\P(t)) = M\P(t)) , (2) 

where the matrix M encodes the transition rates and \P{t)) is the vector of probabilities 
for each configuration. Because M is a stochastic matrix it is guaranteed [32] by the 
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Perron- Frobenius theorem to have eigenvalues satisfying 

A = > Re(Ai) > Re(A 2 ) > . . . . (3) 

The spectrum corresponds to a set of relaxation times t% = — (ReAj)" 1 . The longest 
(non-infinite) relaxation time is r 1; and the associated eigenvalue Ai is the gap which we 
will henceforth denote by the symbol e. At long times the relaxation of any observable 
is expected to decay exponentially with a characteristic timescale T\, a fact we will later 
use to estimate the gap from Monte Carlo simulations. 

The thermodynamic phase boundaries are found by identifying lines along which 
the gap vanishes, indicating the coexistence of two stationary eigenstates (phases). The 
exact calculations of de Gier and Essler [20|l2~T] show that, in the thermodynamic limit 
L — > oo, the gap vanishes along the line a = (3 < ^ that separates the HD and LD 
phases. The gap also vanishes in the entirety of the MC phase, reflecting the generic 
long-range (power- law) correlations seen in this phase [15]. Thus at this level, the static 
and dynamic phase diagrams coincide. 

Where they differ is in the subdivision of the high- and low-density phases^, in 
which the gap remains finite in the limit L — > oo. There is a region, marked LDI' and 
HDI' on Fig. 1(b), within which the gap assumes the asymptotic expression 

2 7T 2 

e(L) = -a-0 + 



(ab)? + l (ab)^-(ab)- 



-L- 2 + 0(L 



-3\ 



(4) 
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1-/3 
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Within the low-density phase (a < /?, a < 
(3 < (3 C where 
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a 



1 — a 



n - 1 

3 



this form of the gap applies for values of 



(6) 



Likewise, in the high-density phase {j5 < a, {3 < |), the region within which the gap is 
given by §4§ is bounded by a < a c where 

(3 



a c (/3) 



P. 



(7) 



In the remainder of the low-density phase, a < > (3 C the gap takes the form 

2 4tt 2 



e(L) = —a — /3 C + - L 

(ab c ) 2 + 1 (ab c ) 2 - (ab c 

Finally, we have by symmetry that when (3 < |,a > a c , 

2 4vr 2 



-L 



-2 



0(L 



-3\ 



e(L) = —a c — (3 + 



[a c by. 



(a c 6)2 - (a c b) 



—L~ + 0(L- 



(9) 



I de Gier and Essler [20H22] refer to these as "massive" phases by analogy with the quantum spin 
chains; we shall only use the terminology associated with the static phase diagram to avoid confusion. 




Figure 3. Exact L —> oo gap (black, dashed) as a function of f3 for a = 0.3, described 
by the HDI' equation ((4]) (red) for j3 < fi c (fl c 0.57, indicated), and by equation flSJ 
(green) for (i > (3 C . 



In these expressions, 

a c = - — — and b c = J^ c . (10) 

« c Pc 

The boundaries between the dynamic subphases — the de Gier-Essler (dGE) lines — are 
shown dotted in Fig. 1(b). Note that the coefficient of the L~ 2 term in this asymptotic 
expansion is discontinuous across the dynamical transition line. 

In order to illustrate the behaviour of the gap along the static and dynamic 
transition lines, we plot it as a function of (3 at some a < \ in Fig. [3J The most 
striking feature is the constancy of the gap above the nontrivial critical point /3 C . Later 
in this work, we will use the constancy of the gap above a critical threshold as an 
empirical means to identify the dynamical transition point. 



2.4- Domain wall theory 

An alternative approach to study the TASEP dynamics is domain wall theory 
(DWT) [33J. Although much simpler than the Bethe ansatz approach of de Gier and 
Essler, it correctly predicts the static phase boundaries and the analytical form of the 
gap in a region of the phase diagram. However, it predicts dynamic subphases that 
are distinct from those found by de Gier and Essler, and that correspond to the static 
subphases. In the numerical study that follows, it will be important to be able to 
distinguish between these two sets of predictions, and so we summarise DWT here. 

The basis of this approach is to assume that collective relaxational dynamics are 
effectively reducible to a single coordinate describing the position of an interface, or 
domain wall. The wall separates a domain of density p~ and current j~ to the left from 
a domain of density p + and current j + = p + (l — p + ) to the right. Each domain is taken 
to possess the steady state characteristics imposed by the boundary reservoir on that 
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Figure 4. The phase diagram obtained from domain wall theory. p + and p~ indicate 
the densities of the left and right domains. In the domain wall theory, static and 
dynamic transitions coincide along the dashed lines. 



particular side of the wall, with p + and p~ therefore as in Fig. HJ The effective theory 
is expected to be exact along the line a = j3 < 1/2 where the exact properties of the 
system such as density profile [15j[T6] and current fluctuations [M] are recovered. 

Let us consider the case a < 1/2 and < 1/2. The motion of the wall is then 
described microscopically as a random walker with left and right hopping rates D~ and 
D + given respectively by imposing mass conservation on the fluxes into and out of the 
wall 125 1 : 



D -_ r _ q(l-a) D+ _ j + _ ^ (n) 

p + — p~ 1 — a — (3 p + — p~ 1 — a — (3' 

For a > (3 the random walk is biased to the left and in the stationary state the domain 
wall is localised at the left boundary and the bulk density is given by p + = 1 — /3. For 
a < (3 the random walk is biased to the right and in the stationary state the domain 
wall is localised at the right boundary and the bulk density is given by p~ = a. Thus 
the first order transition at a = (3 is correctly predicted. 

Moreover, the gap in the resulting spectrum for large L given [25] is given by 

Bdwt = ~D + -D- + 2VD+D- (l - ^ + 0(L- 3 )j . (12) 

Remarkably this expression is identical to (j4j) to order 1/L 2 . Thus in the region 
a < 1/2, [3 < 1/2 DWT correctly predicts the gap. 

When (3 > 1/2 so that p + = 1/2, one can take j + in ( II ip to depend on the size i 
of the right-hand domain. That is, one puts j + — > j + (£) equal to the stationary current 
in a TASEP of size i in the maximal current phase. This implies the large I behaviour 

^4(1 + !' (13) 
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Figure 5. Comparison of the gap for L —> oo from domain wall theory (dashed line) 
with the exact gap (solid line) at a = 0.1 (for which j3 c ~ 0.675). 

This results in a modification of the density profile to an exponential spatial decay 
modulated by a power law with power 3/2, similar to Eq. [1] (LDII). 

In brief, the DWT is remarkably successful, correctly predicting the static phase 
diagram (including subphases), and the exact thermodynamic gap function found by 
de Gier and Essler in the region a < 1/2 and ft < 1/2. It differs in that the dynamic 
subphases are not given by the dGE lines, but the static subphase transition lines. Thus, 
the region within which the gap is constant is different in the two theories, as indicated 
by Fig. El 

3. Evidence for the dynamical transition from Monte Carlo simulations 

We now describe attempts to measure the gap in Monte Carlo simulations of the TASEP, 
with a particular interest in distinguishing between the dGE and DWT predictions. As 
noted previously, the ensemble average of any time-dependent observable 0(t) should 
be dominated at late times by an exponential decay e~' Al ' 4 to its stationary value (O). 
Recall that Ai is the largest nonzero eigenvalue of the matrix M appearing in the master 
equation (j2j), and that all nonstationary eigenvalues have negative real part. In principle, 
therefore, all one needs to do is pick an observable, and measure its late-time decay rate. 
In practice, this is made difficult by the fact that all decay modes may contribute to 
(0(t)), and that by the time the higher modes have decayed away, the residual signal 
(0(t)) — (O) may be extremely small and swamped by the noise. 

Since we are interested in the time- dependence of an observable, we employ a 
continuous-time (Gillespie) algorithm (35] to simulate the TASEP dynamics. More 
precisely, we maintain a list of events (i.e., a particle hopping to the next site, or entering 
or leaving the system) that can take place given the current lattice configuration. 
A particular event i is then chosen with a probability proportional to its rate Ui as 
specified in Section [2j A time variable is then incremented by an amount chosen from 
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an exponential distribution with mean (J2i In this way, a member of the ensemble 

of all continuous-time trajectories of the TASEP dynamics from some prescribed initial 
condition is generated with the appropriate probability. 

3.1. Decay of total occupancy to stationarity 

We consider first the decay of the total number of particles on the lattice, N(t), to its 
stationary value as measured by the function 

R(f) (N(t)) - (AQ 

R{t) ~ (iV(0)> - (N) ■ (14) 
In this expression, angle brackets denote an average over an ensemble of initial conditions 
and over stochastic trajectories. In each simulation run, an initial condition was 
constructed in which each site was independently occupied with a probability p = 1 — /3 
(in the LD phase) or p = a (in the HD phase). In the bulk, these densities then relax 



to the steady-state values displayed on Fig. l(a 



Once R(t) has been sampled over multiple trajectories, the task is to identify a 
time window over which one can fit an exponential decay and therewith estimate a 
gap. The start and end points of this windows are both crucial. If it starts too early, 
then one may expect contributions from subdominant transients (i.e., the decays at rate 
A2, A3, . . .) to systematically skew the estimate of the gap. If it ends too late, noise may 
instead dominate the estimate. 

The noise at the top end we handle by examining the behaviour of local decay rates 

= \nR(t i+i )-\nR(U) 

H+l ~ 

where ti and U + i are successive time points at which R(t) was sampled. At late times, 
one should find ^ — >■ 1. Strong deviations from unity indicate the dominance of 
noise, and we rejected points after which the magnitude of this ratio exceeded a critical 
value. For our data sets, we found that 5 was a suitable choice for this value: see 



Fig. 6(a) for an example. 

The bottom end of the window was chosen by maximising a goodness-of-fit measure 
to a fit of the exponential f(t) = ae~ xt to data points within the window. We adopted 
the adjusted coefficient of determination [36] 

7?2 _ 1 ^(R(u)-f(u)y 

as our goodness-of-fit measure, varying the start of the window until this quantity was 
maximised. In this expression R is the arithmetic mean of R(t) over the set of n times 
ti falling within the window, and k is the number of free parameters in the fit function 
f(t), i.e., k = 2. This goodness-of-fit measure trades off the increased quality of the 
fit obtained by discarding data points against the increasing uncertainty in the fitting 
parameters that comes with the noisier data at the top end of the window. We show 



in Fig. 6(b) an example of how the goodness-of-fit varies with the size of the window. 
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Figure 6. Criteria for choosing the range of data over which to fit a single exponential 



to the decay of the density to stationarity. (a) The ratio of gradients fii+i/fii which 
should be close to unity where a single exponential is a good fit. At large times, this 
ratio becomes noisy; the presence of extreme values (here, a magnitude larger than 5) 
is used to identify the onset of noise. |(b)| The behaviour of 1Z 2 as a function of the 
start of the window. We choose this such that TZ 2 is maximised. In both cases, L = 50, 
a = 0.3 and (3 = 0.64. 



We remark that the relative flatness of TZ 2 above the optimal starting time can be 
ascribed to the fact that only a single decay mode remains in this region. Although this 
procedure is a little cumbersome to apply to each data set, we found it preferable to 
fitting multiple exponentials to the data for R(t), which leads to large uncertainties in 
the measured gap. 

We show results in Fig. in which the gap is plotted as a function of (3 for fixed 
a = 0.3 for a range of system sizes. At the smallest system sizes (L = 10) data 
from Monte Carlo simulation are in agreement with results obtained through exact 
diagonalisation of the matrix M. As the system size is increased beyond the sizes for 
which exact techniques remain tractable, the curves approach those predicted by dGE 
and DWT. However, it is not possible to discern from the largest system simulated, 
L = 150, which of these theories is more appropriate. 

We therefore attempt to extrapolate to the limit L — > oo by fitting a function of 
the form e + g^ -2 + a^L~ 3 + g^L -4 to estimates of the gap at finite system sizes L but 
fixed a and 0. This particular fitting function was found to have the best goodness-of-fit 
(taking into account the varying number of parameters) out of those that were tried, 
and furthermore had the most uniform goodness-of-fit as a function of 0. The results 
of this extrapolation procedure are shown in Fig. El These data appear to exclude the 
functional form for the domain wall theory, although the uncertainties in the estimated 
gaps are such that further evidence is needed in order to confidently rule it out. 

3.2. Decay of other observables 

One possible means to reduce the uncertainty in the measured gap is to try different 
observables 0(t) in the hope that contributions from subdominant decay modes — and 
in particular the second longest-lived mode that decays at rate A 2 — are reduced and 
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Figure 7. The gap as a function of fj at a = 0.3 and system sizes L ranging from 
L = 40 to L = 150 (from bottom to top) obtained from Monte Carlo simulations. The 
L^co dGE (solid line) and DWT (dashed line) predictions are shown for comparison: 
even with the largest system L = 150 simulated, one cannot distinguish between them 
numerically. 
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Figure 8. The gap extrapolated to L = oo from the data shown in Fig. [Jj Again the 
dGE (solid line) and DWT (dashed line) predictions are shown for comparison: here 
the data are potentially more compatible with the former than the latter. 



thereby allow a more accurate estimate of Ai. We considered two candidates. 

The first was the autocorrelation function of the total occupancy, 0(t) = 
N(to)N(t + 1), where t is some fixed time point in the steady state. The analogue of 
the function R(t) is 

(N(t + t)N(t ))-(N) 2 



X ® (N 2 ) - (iV) 2 

where here all averages are taken in the steady state. 



(17) 

Numerically, this quantity is 
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Figure 9. Interactions between first class (circles) and second class (hexagons) 
particles. Second class particles hop can forward into empty sites and exchange places 
with first-class particles behind them. Note that in our simulations, at most one second- 
class particle was present at any instant, and that the behaviour at the boundaries was 
chosen so as not to affect the motion of first-class particles (see text). 



slightly more convenient than N(t), since one can obtain (N(to + t)N(to)) by sampling 
at different to an d t from a single simulation run that has reached the steady state. 

The other observable we investigated was the position of a second class particle 
[37T439] . Like a first-class particle, a second-class particle hops to the right into vacant 
sites. However, it may also exchange places with a particle occupying the site to its left. 
Both of these processes take place at unit rate, and first-class hop to the right at the 
same rate, irrespective of whether the site in front of them is empty or occupied by a 
second-class particle — see Fig. [91 So as not to affect the entry of first-class particles into 
the system, a second-class particle on the left-boundary site is forced to exit if a first- 
class particle attempts to enter. It is reinserted as soon as the left-boundary site becomes 
vacant again. Second-class particles are prevented from exiting at the right boundary. 
In these simulations one can measure the position of the second class particle, x(t), 
starting from an initial condition where each site is occupied by a first-class particle 
with density p (as described above) except for the left boundary site, which is always 
initially occupied by the second-class particle. Here, the analogue of R(t) is 



Gaps obtained from these two functions x{t) an d £,{t) are compared with those 
obtained from R(t) in Fig. [TD] at various system sizes. We see that there are no 
systematic differences between these functions at the system sizes studied, and therefore 
that they are unlikely to provide improved estimates of the gap in the thermodynamic 
limit L — > oo. We therefore conclude that whilst Monte Carlo simulation data is possibly 
more consistent with the de Gier-Essler prediction for the gap than the domain wall 
theory, the effect of the changing gap on the relaxation of a typical observable is very 
small and hard to disentangle from the noise. 

4. Evidence for the dynamical transition from density matrix 
renormalisation 

Given the difficulties in measuring the gap from Monte Carlo simulations, we now 
turn to a fundamentally different numerical approach, namely the density matrix 
renormalisation group (DMRG) [27J that was briefly discussed in the introduction. 
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Figure 10. The gap as a function of as obtained by consideration of the relaxation 
of three different observables. Solid lines show data obtained from the relaxation of the 
density for system sizes L = 10,20,30,50, 150 (bottom to top). Alongside the lower 
three curves are plotted data from the stationary density auto-correlation function (|17p 
(dashed lines) , and alongside the upper two data from the position of a single second- 
class particle (TT5)) (large symbols). All three observables yield consistent results. 



We begin by recapitulating the essential features of the DMRG approach before 
demonstrating that it does indeed show a dynamical transition along the de Gier-Essler 
line. 



4.1. DMRG algorithm for the TASEP gap 

The density matrix renormalisation group (DMRG) procedure builds an approximation 
to the transition matrix M that appears in the master equation (j2J). The basic idea is 
to repeatedly add sites to the system and renormalise the enlarged transition matrix 
so that its dimensionality remains constant as the system size grows. In this way, 
the approximated transition matrix remains sufficiently small that its spectrum (and 
in particular the gap) can be determined using standard numerical diagonalisation 
methods. The success of this procedure relies on finding a reduced basis set at 
each renormalisation step that allows the lowest-lying eigenstates to remain well 
approximated as sites are added. 

Since the TASEP has open boundaries, the extension of the lattice is achieved in 
this instance by adding sites to the middle of the system rather than at one of the 
ends [26]. The system is thus divided into two 'blocks': the left and right halves of the 
lattice. It is extended in each iteration by adding a site at the right-hand end of the 
left block, and the left-hand end of the right block. In order to understand the results 
that come out of the procedure, it is worth explaining in a little more detail how the 



renormalisation is actually achieved — full details are provided in Appendix A 



We outline the procedure in terms of the transformation applied to the left block. 
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Figure 11. Illustration of the DMRG procedure applied to the left block. Before the 
transformation, the state of the block is specified by 1 < p < m and a = 0, 1. After 
adding a new site to the end of the block, the original i sites are renormalised so that 
they are specified by the two numbers 1 < p < m and a = 0, 1. The right block is a 
mirror image of the left, and treated the same way. 



Suppose that at the start of the iteration, the block comprises £ lattice sites. The set 
of states the block can be in is specified by two 'quantum' numbers, p = 1,2 ... , m and 
a = 0, 1. The quantity p indexes the configuration of the first £ — 1 sites, and a the 
occupancy of the rightmost site. The first step of the renormalisation is to construct 
some basis in the 2m-dimensional space spanned by p and a and, crucially, retain only 
m of these basis vectors. This defines a renormalised index p = 1,2, ... ,m for the 
entire £-site block. A new site is then added to the right hand side of the lattice. The 
configuration of this additional site is specified by the renormalised coordinate a. See 
Fig. [TTJ The reason for keeping the rightmost site in the 'physical' (site occupancy) 
basis is that one needs to project the transition matrix for hops from site £ to £ + 1 onto 
the space spanned by p and a. These transition rates are initially only known in this 
physical basis. Thus through this procedure, we obtain a (truncated) representation of 
the transition matrix for an £ + 1-site system in a space of the same dimensionality as 
that used to describe the £-site system. 

The same transformation is applied to the right block, the only difference being that 
it is the leftmost site that is expressed in the physical basis. Interactions between the two 
blocks enter through the construction of the renormalised indices p as we now describe. 
First, a transition matrix M' for the entire system of 2£ sites is constructed by combining 
the matrices for the two halves, and by adding an interaction term that allows a particle 
in the left block to hop to the right block. Again, having the internal two sites specified 
in the physical basis helps here. This transition matrix is then diagonalised. However, 
it is not these eigenstates that are used to perform the renormalisation. Rather, a 
density matrix, which is a symmetric combination of the two longest-lived eigenstates 
of M', is constructed, and eigenvectors of this density matrix are instead used for the 
renormalisation. The idea is that this form of the density matrix allows the stationary 
and longest-lived transient state to be accurately represented in large systems. The 
specific form of the density matrix, and the prescription for obtaining the truncated 
basis set, is given in Appendix A For further details about the principles behind DMRG 
and other applications we refer to |27JllQllIT] , and especially [29|. 
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Figure 12. Stable DMRG data for the gap (red crosses) and the exact thermodynamic 
gap (black line) at a = 0.4 and (3 = 0.75. 



4.2. DMRG results 

We used the DMRG algorithm outlined above to estimate the gap for a given 
combination of a and ft by starting with an exact diagonalisation of the L = 8 system, 
and keeping m — 16 eigenvectors of the density matrix in each renormalisation step. In 
principle, one ought to be able to access arbitrarily large system sizes with this method. 
In practice — and as was also noted in [26] — the algorithm eventually goes unstable, 
which is typically manifested through the gap acquiring an imaginary part. We simply 
ignore data for system sizes where the instability is judged to have kicked in. 

Although we can access larger system sizes with the DMRG approach than was 
possible with Monte Carlo (e.g., L up to about 250, as shown in Fig. fl2l) . it is still 
necessary to extrapolate to the thermodynamic limit, L — > oo. We follow a similar 
procedure to that described in Section 13.11 That is, we specify a finite-size fitting 
function of the form f(L) = e + g^ -2 , and adjust the smallest value of L used for the 
fit. Again, we use V? as a goodness-of-fit measure, i.e., Eq. ({16]) with t replaced by L. 



We show in Fig. 13(a) how the goodness of fit varies with the minimum value of L used 



in the fit; choosing the optimal (largest) value yields an estimate of the gap that we 



observe from Fig. 13(b) that appears more consistent with the de Gier-Essler prediction 
than domain wall theory. 

We show the DMRG estimate of the gap obtained in this way as a function of ft 
in Fig. [TH The agreement with the analytical prediction common to dGE and DWT 
below ft < \ is very good. For larger ft the data are scattered around an apparently 
constant value, indicating the presence of inaccuracies in the DMRG algorithm and/or 
the extrapolation of L — > oo. If we regard these data as independent measurements of 
the same value, we can obtain an estimate of the uncertainty in the constant value of 
the gap above some critical point by simply calculating the standard deviation. We find 
the DMRG gap to approach the constant value e = —0.00125(9) when a = 0.4. For 
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Figure 13. (a) Coefficient of determination R and |(b)| extrapolated gap, as a 
function of I/ m i n for DMRG data corresponding to a — 0.4, f3 = 0.75. The optimal 
choice of L m i n lies just below 200, for which the extrapolated gap very closely matches 
the de Gier-Essler result and clearly differs from the DWT gap. 
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Figure 14. DMRG estimates of the gap (red crosses), exact thermodynamic gap 
(black curve) and DWT gap (green curve), all as a function of (3 at a = 0.4 (for which 
p c « 0.53). 



comparison, de Gier and Essler predict a constant value for the gap of —0.00121 . . . above 
the dynamical transition, while domain wall theory predicts —0.00102 . . .. The DMRG 
measurement is then clearly consistent with the de Gier-Essler prediction, whilst the 
difference from the domain wall theory value is significant. Taking this result together 
with the Monte Carlo data we conclude that a dynamical transition does indeed occur 
at the dGE line, and not where DWT predicts. 

5. Discussion of the distinction between static and dynamic subphases 

In this paper we have presented numerical evidence, the most convincing of which came 
from DMRG calculations, that a dynamical transition occurs in the TASEP at the 
dGE line rather than at the subphase boundary. However, this finding in turn raises 
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a number of open questions. For example, does the dGE line — and in particular its 
departure from the static subphase boundary — have any physical significance? How is 
the static subphase boundary manifested in the eigenvalue spectrum? We discuss the 
distinction between the two phase diagrams first in general theoretical terms, and then 
with reference to different treatments of the TASEP. 

5.1. General theory of nonequilibrium phase transitions 

One way to gain a general understanding of nonequilibrium phase transitions is through 
the analogue of partition function for a system of size L, Zl [31,32j. It has been shown 
that for any system governed by Markovian dynamics, it may be written as 



where (— Xj) are the eigenvalues of the Markov transition matrix which defines the 
dynamics [3"2"jr4*2] . For a finite, irreducible configuration space there is a unique stationary 
state corresponding to the largest eigenvalue A = 0. By the Perron-Frobenious 
theorem [43], Ao is non-degenerate, therefore the gap is given by the second largest 
eigenvalue Ai and the largest relaxation time is T\ = — 

From our knowledge of equilibrium phase transitions, we expect a static phase 
transition to occur when lim^oo \ In Z L exhibits nonanalyticities. We see from ffT9"j) 
that to obtain a static phase transition we require a major restructuring of the eigenvalue 
spectrum. To be general and explicit let us consider the nonequilibrium analogue of the 
free energy density, 



Consider an arbitrary eigenvalue Aj in ( l2"Uj) . possessing a nonanalyticity at a critical value 
of some parameter. Unless a significant number of other eigenvalues in the spectrum 
converge onto the same nonanalyticity in the thermodynamic limit, the effect of the Aj 
nonanalyticity is lost as L — > oo and (1201) and hence hi remain analytic. 

On the other hand a dynamical phase transition, as defined above, is only concerned 
with the eigenvalue Ai. Therefore a dynamical phase transition does not necessarily 
coincide with a static phase transition. One very simple example would be if the two 
leading subdominant eigenvalues Ai and A2 crossed at some values of a control parameter. 
Then the gap would behave in a nonanalytic way but Z would be analytic. 

Let us consider more explicitly the TASEP. From the exact solution of the TASEP 
[4"]IT5] the expression for Z L , as defined by ( IT9|) . is 




(19) 




(20) 




(21) 



where 




(22) 
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and for L large 



X < 



1 



x(l - 2x) 

[x(l - X)] L + X " 2 



'ttL 

xA L 



x 



X > 



I • ( 23 ) 



0F(2x - 1) 2 L 3 / 2 

Now consider, for example, the case a < 1/2. As we increase /3 a first-order static 
phase transition occurs on the line ft = a where the dominant contribution to (j2ip 
changes from $l(/3) t° < ^ ) l(q ; )- Thus at the phase transition h L defined in f )20|) is non- 
analytic. The subphase boundary is at ft = 1/2 where the asymptotic behaviour of 
the subdominant contribution, changes. As this is a subdominant contribution 

to InZi, hi remains analytic at the sub-phase boundary. At the dGE line there is no 
noticeable change in the asymptotic behaviour of and is analytic. 

5.2. Static vs dynamic transitions within various treatments of the TASEP 

We illustrate further the distinctions between the static and dynamics phase diagrams 
by collecting exact results for the TASEP together with the predictions of various level 
of approximations such as mean field theory and domain wall theory. 

Burgers Equation (Mean Field theory) We begin by considering a mean field 
description of the system that is given by the Burgers equation, 

d tP (x,t) = -(l-2p)d x p+^d 2 x p, (24) 

where p(x, t) is a density field. Although the Burgers equation has been studied 
extensively in the literature, we are not aware of any treatment of the case with 
prescribed boundary reservoir densities. In the appendix we give the solution of (1241) 
subject to the boundary conditions p(0) = a and p(L) = 1 — (3, and arbitrary initial 
condition. 

The results for the Burgers gap eb may be summarised as follows (restricting 
ourselves to the case a < 1/2): 

For /? < 1/2 e B = - 2\a(l - a) - p(l - p)\ (25) 
For > 1/2 e B = - (1 ~ 2 2a)2 (26) 

The behaviour of the gap is plotted in Fig. [15j Note that a dynamical transition occurs 
at (3 — 1/2. Interestingly, it turns out that there is no change in the density profile at 
this value so, in fact, the mean field theory predicts a dynamical transition instead of a 
subphase boundary at ft = 1/2. Thus even at the level of mean-field theory, the TASEP 
exhibits the phenomenon of a dynamic transition that is not accompanied by a static 
transition. 
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Figure 15. Comparison of the gap for L — > oo from the spectrum of the viscous 
Burgers equation (l24l) to the exact thermodynamic gap, for a = 0.1 and f3 c « 0.675. 

Domain wall theory The predictions of domain wall theory were reviewed in 
Section 12.41 As we pointed out there, DWT gives a subphase boundary at (3 = 1/2 
along with a dynamical transition: the distinction between the two phase diagrams is 
lost. However, in the subphases LDI and HDI, the DWT gives the correct dynamical 
gap to order 1/L 2 . 

Exact results As was discussed in Section |2} the exact results for the stationary state 
[T5|[T6] and the spectrum [20~|l2Tj reveal behaviour that is different from both mean-field 
and domain wall theory. In common with the latter, there is a static subphase boundary 
at = |. Like the mean field theory, the location of the dynamic transition occurs at a 
distinct location to the static phase boundary. However, this location a c {j5) is nontrivial 
and is not predicted by either of the simpler theories. Moreover, the behaviour of the gap 
at the static transition point f3 = a differs between the mean-field and exact theories: 
in the former, the gap has a cusp, whereas in the latter it varies analytically. 

6. Conclusion 

In this work we have provided numerical evidence to confirm the existence of the 
dynamical transition line predicted by de Gier and Essler. Ultimately the evidence 
came from the approximate, but quantitatively reliable, DMRG technique. 

To conclude, we now return to the main open question resulting from this work, 
namely that of the physical significance of the dynamical transition that takes place 
along the de Gier-Essler line. We remark again that, as was seen in Section [31 
direct measurement of the gap in Monte Carlo simulations was very difficult. The 
nonanalyticity in the leading eigenvalue was barely detectable in the stochastic 
simulation data. We were therefore unable to identify, for example, whether the system 
relaxed to stationarity in a fundamentally different way either side of the transition line. 
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Given the success of the domain wall theory in predicting the gap in part of the 
phase diagram, it is tempting to believe that with some refinement it may be able to 
explain the dynamical transition and furthermore reveal the physics associated with it. 
For example, it has been noted in [21], that DWT may be modified to give the correct 
gap throughout the LD1' region by simply imposing that the right domain density be 
p + = 1 — (3 C in this region. However DWT then no longer predicts a static subphase 
boundary at (3 = 1/2 and also the transition line a c ((3) is not predicted. 

It may be that some knowledge of the form of the eigenvectors associated to the 
low lying states will allow an understanding of what occurs at the dynamical transition. 
It would be of interest to extend the DMRG studies to investigate the form of the 
eigenfunctions. Also recent analytic work by Simon jH], in which eigenvectors are 
constructed for the partially asymmetric exclusion process through the co-ordinate Bethe 
ansatz, may help to shed light on the matter. 
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Appendix A. DMRG algorithm for TASEP 



In this appendix, we provide a step-by-step description of the DMRG algorithm used 
to obtain the results presented in Section HI 

In what follows we make use of three elementary transition matrices for the TASEP. 
With the ordering (_ , • ) of single-site configurations we define 



hi 



-a 
a 



and h r 



(Al) 



(^) 

for the entry and exit of particles at the left and right boundary sites. In the bulk, and 
with the ordering ( ,_•>•_,••) of two-site configurations, we define 

/ \ 

0010 /A o^ 

h b= n n in " ( A - 2 ) 



( 








\ 








1 











-1 





V o 








o J 



The DMRG algorithm now proceeds as follows. 

0. Initialise block transition matrices We partition the system into two blocks, each 
of size i. The transition matrix for the left block can be written as 

e-i 



M. 



hr®J^ 



i-i 



E^ 1 

k=l 



® fr 6 ® X 5 



l-k-l 



(A.3) 
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where X is the 2x2 identity matrix. Likewise, for the right block one has 

e-i 

Mj? = ]T X®^ 1 ® h b <g> X®^" 1 + X^" 1 ® /i r . (A.4) 
fe=i 

1. Diagonalise the transition matrix for the entire system The transition matrix for 
the entire system of size 21, is given by 

M (2Q = M W g j®« + x ^-i g ^ g jf-i + x i g M W (A-5) 

Now we diagonalise M^ 2e \ first using the Arnoldi algorithm |35] (as implemented 
in the ARPACK package jl6] and accessed using Mathematica) to find the left 
and right ground state eigenvectors (<f>o\, \if>o)- We also want to compute the gap 
and associated eigenvectors \ipi) and Following [26], we find this can be done 
more accurately by 'shifting' the transition matrix through M' = + A|^o)(0o| 
sufficiently far that the gap is now the largest eigenvector of M' . 

2. Form reduced density matrices The reduced density matrices are defined as 

PL =\Tr R ^){^\ + |0 o )(0o| + IViX^I + \<h)(<h\) (A.6) 

p R = ^Tr L (\^ )(H + |0o)(0o| + MM + |0i)(0i|) • (A.7) 

Here Tr^ and Ttr indicate a trace over the degrees of freedom in the left and right 
blocks respectively. 

To be clear, let and {\j) R } be basis sets for the left and right blocks 

respectively. A state vector \ip) across both blocks can then be expanded as 

H) = H^l\j)r. (A.8) 
The trace operation Tr£,(|^>)(^|) is then defined as 

Ttdmm = E fe<w) \j)r(j'\r • (A.9) 

jj> \ i J 

Note that this is an operator acting only on the right block. The corresponding 
operation Tr R is defined analogously. 

3. Diagonalise the density matrices and form a truncated basis set We first find all 
eigenvalues and eigenvectors of the symmetric matrices pl and pr using a dense 
matrix diagonalisation routine from the LAPACK package [17] , again accessed using 
Mathematica. We then form the matrices Ol and Or that have as their columns 
the m eigenvectors with largest eigenvalues of pi and pr respectively. Then we 
construct projectors 

P L = L ®X and P R =X®0 R (A.10) 

that will be used in the renormalisation step below. 

4. Enlarge system The left block is enlarged by adding a site at its right end. The 
transition matrix for this enlarged block is 

M L e+1) = Mf ® X m + X^" 1 <g> h b . (A. 11) 
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Likewise, the right block is enlarged by adding a site at its left: 

M { 1 +1) = h b ® X"- 1 + Z® 2 ® Mf . (A. 12) 

5. Renormalise the transition matrices for the enlarged system The renormalisation 
is performed by projection the enlarged transition matrices onto the subset of 
density matrix eigenstates that was retained in step 3. Formally, this is achieved 
by computing 

Mt 1] = PlMt 1] P L and H { * +1) = P T R M { £ +1) Pr • (A.13) 

Note that after this transformation, the first i sites of the system are no longer 
represented in the 'physical' basis (i.e., particle-hole configurations). However, the 
form of the projectors (1A.10I) ensures that the last site of the block is represented 
in the physical basis, which is essential for the expression (I A. 5 1) to be valid at each 
stage of the renormalisation. 

6. Return to step 1, putting £ — > t + 1. 

Appendix B. Solution of the Burgers equation with fixed boundary 
densities 

Here we solve the Burgers equation (i.e., the continuum limit of the TASEP within a 
mean-field theory) 

d tP {x,t) = -{l-2p)d x p+ l -dlp. (B.l) 

subject to the boundary conditions p(0) = a and p(L) = 1-/3. We use the Cole-Hopf 
transformation |48j, which involves the change of variable 



9 , , . 
1 + —mu(x,t) 
ox 



(B.2) 



: u(x,t) = -— u(x,t) . (B.3) 



Substituting into (IB. II) . and integrating with respect to x eventually yields 

d , > Id 2 

Wt u ^ t) = 2Tx 2 

The boundary conditions on p(x, t) transform to the conditions 

u'{0,t) = (2a-l)u(0,t) (B.4) 

u'(L, t) = (1 - 2fi)u(L, t) . (B.5) 
We proceed by finding the eigenfunctions (f> n (x) that satisfy 
1 d 2 

;4>n(x) = X n 4>n{x) (B.6) 



2 da; 2 

along with the above boundary conditions at x = and x = L. Thus given an initial 
condition p(x, 0) the function u(x,t) will be given by 

u(x,t) = ^c n e A "V n (x) (B.7) 

n>0 
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where the coefficients c n will depend on the initial condition. Inverting the Cole-Hopf 
transformation gives 

3 (A„-Ao)t . 



p(x,t) 



0q(x) Dn>0 C « e 

0o(x) 



Wo (*) J 



E„> Cne( A "- A »)Vn(^) 



(B.8) 



The final term, involving the ratio of sums, vanishes as t — > oo, as all A n with n > 
exceed A . Thus the stationary profile is 



p(x) 



1 + 



00^) 



(B.9) 



We can define the Burgers gap, Eb via the asymptotic rate of the exponential decay 
of the density profile to its stationary form, viz, 



(B.10) 



e B = Jim — In [p[x, t) - p{x)\ = A x - A . 

t-*oa dt 

The eigenfunctions given by ( IB. 61) themselves take the form 

= Ae 7X + Be~^ x (B.ll) 

whence A = \^ 2 . The allowed values of 7 are quantised due to the boundary conditions 
flBTil) and f lR5l) . We first look for solut ions with A > 0, i.e., real 7. Imposing the 
boundary conditions leads to the pair of equations 

-f(A-B) =(2a-l)(A + B) (B.12) 

7 (Ae 7i - Be~^ L ) = (1 - 2(3){Ae^ x + Be~^ x ) . (B.13) 

These equations are consistent only if 

2(l-a-/3) 7 



tanh(7L) 



(B.14) 



(l-2a)(l-2/3)+ 7 2 

In the limit L — > 00, the left-hand side approaches ±1. The resulting equation has a 
solution 7 = 1 — 2ctifa:<! and a solution 7 = 1 — 2/3 if (3 < |. 

Solutions with a negative A have purely imaginary 7 = ifc (k real), and by following 
through the same procedure, one finds that the allowed values of k are the roots of 

2(o + /3- l)k 



tan(kL) 



(B.15) 



k 2 - (2o- l)(2/3- 1) ' 

In the large-L limit, we find that k ~ 1/L, and hence the associated eigenvalues vanish 
as 1/L 2 in the thermodynamic limit. 

We can now state the L — > 00 behaviour of the Burgers gap e as a and /3 are varied. 
When both a and (3 are smaller than |, there are two isolated positive eigenvalues 
A = ^ 2a ~^ , ^ 2/3 ~ 1 - ) . In this region, then, the gap behaves as 



|(2a-l) 2 - (2/3-1) 



e b = -^ — = -2|o(l-a)-/3(l-/3)| . (B.16) 

Like the exact gap 01]), this vanishes along the HD-LD coexistence line o = /3 < |. 
However, unlike the exact gap, the mean-field gap has a nonanalyticity along this line. 
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If a or is increased above ~, one of the two solutions with positive A ceases to 
exist, and the second-largest eigenvalue Ai = in the thermodynamic limit. Hence then, 
the size of the gap is simply equal to the value of the largest eigenvalue. That is, within 
the low-density phase a < ~, /? > ~, 

EB = (B.17) 
and within the high-density phase a > ^, < \ 

£B = ■ (B.18) 



2' 



The Burgers gap thus exhibits nonanalyticities along the coexistence line a = j3 < 
and along the lines a < = ^ and a = < |. Thus the dynamic phase diagram 
for the Burgers equation appears the same as the static phase diagram, Fig. 1(a) , except 
that the subphase boundaries have become dynamical transition lines. 
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